#!/usr/bin/env perl

=head1 NAME

lcs.pl

=head1 SYNOPSIS

lcs.pl pdbID [pdbID] [-r] [-k] [-a] [-c] [--pathCandidates] [--pathCofm] [--pathResult] [--log]

=head1 DESCRIPTION

Largest common subcomplex (lcs) searches for shared subcomplexes of a given protein complex with other complexes.
The first parameter specifies the PDB-ID of the first complex. If it is followed by a second PDB-ID, only these
two complexes are compared. Otherwise, the first complex is compared to all other complexes. The algorithm works 
by first identifying small shared subcomplexes (of only two domains) and then merging these into larger subcomplexes.

=head1 OPTIONS

=head2 -r maximal RMSD that must not be violated by any subcomplex (default value 5)
=head2 -k return only best 'k' results based on RMSD (default value inf = return all results) 
=head2 -a approximationThreshold triggers the use of an approximative heuristic to limit the processing time. 
Specifically, 'a' limits the number of subcomplexes that are passed to the next merging stage. (default value inf = do not use approximative heuristic)  
=head3 -c compare the query only to the 'c' best suited candidate complexes (default value inf = compare to all suited candidate complexes)
=head4 --pathCandidates path to the 'candidate' file that list the best suited candidate complexes for any query (default "./candidates")
=head5 --pathCofm path to cofm-files that have the coordinates and classifications of all domains (default "./scop_only/")
=head6 --pathResult path to which the result files are written (default "./")
=head7 --log a log file is created (under 'pathResult'), if this option is enabled

=head1 REQUIRES

sbglib, PDL

=head1 AUTHOR

Philipp Harth <phil.harth@gmail.com>

=cut

use Getopt::Long;
use Pod::Usage;
use Log::Log4perl qw(:easy);

use QueryProcessor;
use Scorer;
use ResultWriter;
use Util;


#set initial values for parameters and options
my $pdbId1 = ""; #required parameter
my $pdbId2 = "all"; #optional parameter
my $maxRmsd = 5; #optional
my $bestK = inf; #optional
my $approximationThreshold = inf; #optional
my $bestCandidates = inf; #optional
my $pathCandidates = "./candidates"; #optional
my $pathCofm = "./scop_only/"; #optional
my $pathResult = "./"; #optional
my $log = 0; #optional

#check and load options
GetOptions("r=f" => \$maxRmsd,
		   "k=i" => \$bestK,
		   "a=i" => \$approximationThreshold,
		   "c=i" => \$bestCandidates,
		   "pathCandidates=s" => \$pathCandidates,
		   "pathCofm=s" => \$pathCofm,
		   "pathResult=s" => \$pathResult,
		   "log" => \$log );

#check and load parameter(s)
if(scalar(@ARGV) == 0){
	pod2usage({-message => "No pdbID specified"}); 
}
if(scalar(@ARGV) >= 3){
	pod2usage({-message => "More than 2 pdbIDs specified (maybe you forgot '-' before an option ?)"}); 
}
#read pdbid 1
if(length($ARGV[0]) != 4){
	pod2usage({-message => "pdbID '$ARGV[0]' is invalid"});
}
$pdbId1 = lc $ARGV[0];
#read optional pdbid 2
if(scalar(@ARGV) == 2){
	if(length($ARGV[1]) != 4){
		pod2usage({-message => "pdbID '$ARGV[1]' is invalid"});
	}
	$pdbId2 = lc $ARGV[1];
}

#init logger if '--log' option set
if($log){
	my $logfilepath = join '',">>", $pathResult, $pdbId1,"_log"; 
	Log::Log4perl::easy_init({	
		level => $INFO,
		file => $logfilepath	
	});
}

my $logger = Log::Log4perl::get_logger();
$logger->info("$pdbId1 vs. $pdbId2: r=$maxRmsd, k=$bestK, a=$approximationThreshold, c=$bestCandidates");
$logger->info("filepaths: pathCandidates=$pathCandidates, pathCofm=$pathCofm, pathResult=$pathResult");

#initialize modules
my $scorer = Scorer->new(maxRMSD=>$maxRmsd);
my $queryProcessor = QueryProcessor->new(scorer => $scorer, pathCofm => $pathCofm);
my $resultWriter = ResultWriter->new(filepath => $pathResult);


if($pdbId2 eq "all"){ #query vs. all

	my @candidates = Util::getCandidates($pathCandidates, $pdbId1);
	#repeat for the $bestCandidates first candidates
	for(my $i=0; $i<scalar(@candidates); $i++){
		INFO("$pdbId1 vs. ",$candidates[$i]);
		my @subcomplexes = ();
		@subcomplexes = executeQuery($pdbId1,$candidates[$i],$bestK);
		$resultWriter->write(\@subcomplexes, $pdbId1, $maxRmsd, $bestK, $approximationThreshold, $bestCandidates);
			
		last if $i+1 >= $bestCandidates; #option (default value inf = compare to all candidates)
	}

} else { #one on one matching

	my @subcomplexes = ();
	@subcomplexes = executeQuery($pdbId1,$pdbId2,$bestK);
	$resultWriter->write(\@subcomplexes, $pdbId1, $maxRmsd, $bestK, $approximationThreshold, $bestCandidates);

}

=head2 executeQuery

 Function: implements the query loop
 Example : 
 Returns : list of subcomplexes
 Args    : PDB-ID1, PDB-ID2
 
=cut
sub executeQuery{
	my ($pdbId1, $pdbId2) = @_;	
	#reset query processor, this function may 
	$queryProcessor->reset();
	my $l=1;
	my @results = ();
	#query with mergeOps set to inf. Here 'approximationThreshold' is used, rather than limiting the number of merge operations.
	#This means: do not interrupt merging, until 'query' signals that merging is complete
	#However, as 'approximationThreshold' is other than inf, 'query' may use an approximative heuristic to finish merging
	while($l){
		my ($loop, @intermediateResults) = $queryProcessor->query($pdbId1, $pdbId2, $bestK, inf, $approximationThreshold);
		$l = $loop;
		push @results, @intermediateResults;
	}
	return @results;
}